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ABSTRACT 

Differential  cuadrature  is  a  useful  numerical  technique 
for  solving  non-linear  partial  differential  equations.   It 
involves  approximating  the  partial  derivatives  by  a  linear 
combination  of  functional  values  and,  therefore,  provides 
an  easy  method  of  transformation  of  partial  differential 
equations  into  a  set  of  ordinary  differential  equations. 
The  technique  is  employed  for  solving  boundary  value  problems 
which  can  be  represented  by  partial  differential  equations. 

Most  other  methods  like  the  finite-difference  method 
involve  approximation  in  terms  of  functional  differences 
instead  of  functional  values  and  therefore,  require  functional 
evaluation  at  a  large  number  of  points  for  satisfactory  re- 
sults.  It  is  in  this  respect  that  differential  quadrature 
has  its  major  advantages  over  other  methods  in  terms  of  both, 
the  computer  storage  and  computational  time.   However,  the 
success  of  the  method  depends  largely  upon  the  method  of 
evaluation  of  weighting  coefficients.   Three  methods  are  con- 
sidered in  this  respect  viz. classical  quadrature  analogy, 
Legendre  polynomial  approach  and  spline  approximation. 

Differential  quadrature  is  applied  to  solving  several 
models  in  engineering  with  both  fixed  and  moving  boundary 
conditions.   A  moving  boundary  condition  is  specified  at  a 
■ooint  which  itself  varies  as  a  function  of  time.   Differential 
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quadrature  is  used  to  solve  the  isothermal  reactor  model  as 
well  as  the  adiabatic  reactor  model.   A  lot  of  computer  memory 
and  computation  time  are  saved  "by  using  this  technique. 


CHAPT2H  I 
DIPPER3R2IAI  C.UADRATUR3  AND  SPLIK3  APPROXIMATION 

-  A  LITERATURE  SURVEY 

1.1  Introduction 

Fartial  differential  equations  are  frequently  encountered 
in  the  fields  of  engineering  and  science.  However,  the  numeri- 
cal solution  of  time-dependent  non-linear  partial  differential 
equations  has  been  a  complicated  and  highly  problem  dependent 
process.   In  general,  the  solution  of  slightly  different  type 
of  partial  differential  equations  may  require  separate  and 
completely  different  computer  programs.   Thus,  effective 
numerical  technique  in  this  respect  can  be  very  beneficial. 

The  conventional  numerical  techniques  such  as  the  method 
of  lines,  finite  difference  method,  etc.  require  the  function 
value  to  be  evaluated  at  a  large  number  of  points  to  obtain 
satisfactory  results.   This  requires  a  lot  of  computer  stor- 
age and  computer  time  and  thus  an  increased  cost  and  effort. 

Quadrature  techniques  like  trapezoidal  rule,  Simpson's 
rule  have  been  used  since  early  times  to  estimate  the  area 
under  curves.   A  more  common  technique  is  that  of  Guassian 
quadrature  which  can  •orovide  a  good  approximation  for  inte- 
grals.  This  chanter  deals  with  a  recently  developed  numerical 
technique  loiown  as  differential  quadrature  which  is  very 
much  similar  in  nrinci-cle  to  that  of  Guassian  cuadrature  and 
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can  be  effectively  use  a  for  polynomial  approximation  of 
partial  differentials.   The  follovring  chapters  deal  with  some 
of  the  applications  of  this  method  in  the  field  of  engineer- 
ing. 

Consider  the  type  of  approximation 

Lu  =  y     a.  u.  (1) 

i 

which  is  l:nor,n  as  quadrature  approximation  if  L  is  an  inte- 
gral o-cerator.   In  analogy,  the  same  approximation  v,-as  named 
as  differential  quadrature  "by  Bellman  (1)  when  L  is  an  inte- 
gral operator.  Unlilie  Guassian  quadrature,  differential 
quadrature  method  is  still  at  an  early  stage  of  development. 
The  advantages  of  the  method  will  be  obvious  after  the  dis- 
cussion of  the  follovdng  chapters. 

Differential  cuadrature  method  is  different  from  most 
conventional  methods  in  the  sense  that  the  interpolation  is 
expressed  in  terms  of  the  values  of  the  function  instead  of 
the  differences  of  the  function.   However,  the  approximation 
of  derivatives  by  differences  provide  the  basis  of  many  me- 
thods of  solving  differential  equations.   Interpolation  using 
functional  values  at  certain  selected  points  has  been  dropped 
in  past  due  to  hazards  of  roundoff  errors,  but  Bellman  and 
his  coworkers  has  shown  that  the  difficulty  can  be  overcome 
by  first  smoothening  the  data  and  then  differentiating.   They 
have  successfully  applied  the  method  to  the  solution  of  fluid 
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flow  equations,  the  Kodgkin-Haxley  model  (2)  and  other  models. 

1.2  Differential  Quadrature 

Consider  a  linear  or  non-linear  first  order  partial  dif- 
ferential equation  of  the  form 

Ut  =  f(t,  x,  U,  U„)  (2) 

where  x  lies  in  the  finite  interval  (a,b)  and  the  boundary 
conditions  of  the  problem  are  in  any  given  form.   The  initial 
condition  is  assumed  to  "be  knovoi  and  is  of  the  form 

U(o,x)  =  U°(x)  (2) 

Assuming  the  function  U  to  he  sufficiently  smooth  in 
the  interval  (a,b)  we  can  write  the  following  approximate 

solution  (3) 

U„  (t,  x.)  =  I   a..  U  (t,x.)  ,    =  1,  2,  .  .  .  ;N 

(A) 

where,  N  is  the  number  of  mesh  points  selected  and  a^.  is  the 
matrix  of  weighting  coefficients  of  order  NX!!.   N  is  also 
known  as  the  order  of  differential  quadrature  method.   Sub- 
stituting equation  (4)  into  equation  (2),  we  get, 

U(t,xi)  =  f(t,xif  U(t,x±),  £   a±.  UCt.^W  , 

3=1 

i  =  1,  2,  .  .  .  ;1"I      (5) 
which  is  a  set  of  IT  ordinary  differential  equations  with 
initial  conditions 

U(0,X±)  =  U°  U±),    i  =  1,  2,  .  .  .  ;IT     (6) 
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Thus,  the  solution  of  equation  (2)  can  be  obtained  by 
solving  the  set  of  ordinary  differential  equations  (5)  v.lth 
initial  condition  (5).   The  boundary  conditions  for  eon. (2) 
can  be  reduced  to  a  set  of  algebric  equations  by  using  ap- 
proximation (2)  and  the  result  can  be  used  in  order  to  elimi- 
nate two  of  the  IT  variables,  U(t,x.),  i=l,  2,  ,11  in 

equations  (5).   Thus,  the  solution  is  obtained  by  solving 
the  (IT-2)  ordinary  differential  equations  and  the  remaining 
two  functional  values  can,  therefore,  be  evaluated  using 
the  boundary  conditions  of  the  problem. 

For  approximation  of  higher  order  methods,  the  same  idea 
can  be  extended.   Considering  equation  (2)  as  a  linear  trans- 
formation of  U,  we  can  v.*rite 

Ux  =  AU  (7) 

Then,  the  second  order  derivatives  can  be  approximated  as 

U.„._  =  (U,,)x  =  A  (AU)  =  A2U  (8) 

The  higher  order  derivatives  can  similarly  be  found  by 
iterating  the  linear  transformation  A.   Writing  equation  (8) 
in  the  sense  of  equation  (4),  we  have, 

U_(t,::,  )  =  X    i       -ik  *}  .  U   (t), 
k=l   3=1 

•i  -  1         0  V  ( Q ) 

Similarly, 

N   II  IT 
U.__(t,r..)  =Z  Z  Z        a.,  ar  a,   U,(t), 
x   i=l  k=l  j=l   —  —  "d   J 

i  =  1,  2,  ,1T        (10) 
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anc.  so  on. 

Therefore,  it  can  "be  seen  that  success  of  differential 
quadrature  method  depends  largely  upon  the  values  of  weighting 
coefficients  a.  . *  s.   The  method  which  have  been  used  in  the 
•oast  by  Bellman  and  his  coworkers  (4,5,6)  are  discussed  in 
the  next  section. 

1.3  Determination  of  weighting  Coefficients  a^.'s 

1.3.1  3y  analog5/  with  the  classical  quadrature  case 

The  coefficients  a.  -'s  in  the  approximation, 

U„(xs  )  =  -y     a.  •  U(x.  )j   i  =  1,  -,  .  .  .  .  ,-• 
3=1  (11) 

can  be  easily  determined  by  analogy  with  the  classical  quadra- 
ture case  which  demands  that  equation  (2)  be  exact  for  all 
polynomials  of  degree  less  than  or  equal  to  (H-l).   Consider- 
ing the  test  function, 

g1:(x)  =  x1:,     k  -  1,  2,  .  .  .  .  ,N    (12) 

For  arbitrary  points  :•:.,  ,  i=l,  2,  .  .  .  .  ,N  this  leads  to 
a  set  of  linear  algebric  equations 

£  a,  .  x,-1  =  (k-l)  x  *   X  ,   k  =  1,  2,  .  .  .  .  ,1! 

3=1   D   3 

i  -  1       2  N 

j-  —  j- ,  ^ ,  .  •  .  .  ,i.> 

(13) 

Thus,  choosing  I!  and  %..  ,    i  =  1,  2,  .  .  .  .  ,N  the  values 

0^  the  reis:htin£  coefficients  a.  .'s  can  be  unicuely  determined, 

"  13 
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1.3.2  By  analogy  with  Lagrange's  interpolation  formula 

Instead  of  solving  a  set  of  algebric  equations  in  the 
previous  case,  the  coefficients  a. .*s  can  also  he  determined 
by  pronerly  selecting  x.'s  i.e.  if  we  consider  :■:  to  be  the 
root  of  shifted  lagendre  polynomial  of  degree  N,  ?«*  (xj 
with  orthogonality  range  of  (0,1),  the  legendre  polynomial 
?„*  (::)  being  defined  as 

?  *  (::)  =  ?TT  (l-2x)  (14) 

3y  analog?/  v."ith  lagrange's  interpolation  fomula  the 
trial  function  is  taken  of  the  form, 

gk(x)  =  PK*  (::)/  [(::-::.)  PH*  •  (.j]    (15) 

Since  P.T*  (:•:)  is  a  polynomial  of  degree  IT,  g,  (x)  will 
be  a  polynomial  of  degree  (K-l)  such  that 

Sk(xj)  =  bih  »   h  =  1,  2,  .  .  .  .  ,K 

j  =  1,  2,  .  .  .  .  ,1T    (16) 

Assuming  that  the  equation, 

Ux(::i}  =X  aij  U(xj)j  i  =  1»  2'  '  *  '  *  ,IT 
^=i  (17) 

is   exact   for  U (:•:)  =  gv (:■:) ,    we   have 
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lZ2l (18) 

("±   -  x,J   ?N*(x,J 


For  i=k  we  use  the  L* Hospital  rule  and  the  fact  that 
?  *(x)  satifies  the  differential  equation 


x(l-x2  )v  *»    [x)  +  (l+2x)?  *«(x)  +  N(N-l)Pw*(s)  =  C 

(19) 
whi  ch  give  s , 

a^  =  d-2r1:)/  [2::,.( ::,.-!;]  (20) 

Thus  the  weighting  coefficients  a.  .*s  can  be  calculated 
by  equation  (18)  and  (20). 

A  better  approach  to  the  determination  of  weighting 
coefficients  is  given  by  snoothening  the  data.   Spline  ap- 
proximation is  one  of  the  best  methods  in  this  respect.   It 
is  discussed  in  the  following  section. 

1.4  Spline  Approximation 

1.4.1  Introduction 

A  srline  in  its  simplest  form  can  be  understood  from 
analogy  with  a  draftsman  srline  which  is  a  thin  strip  of  wooc 
or  some  other  material  used  for  fitting  smooth  curves  passing 
through  certain  -points.   These  splines  are  anchored  in  place 
by  applying  lead  weights  called  "ducks"  at  points  along  the 
srsline  and  the  snline  can  thus  be  made  to  pass  through  cer- 
tain -ooints  by  adjusting  the  position  of  these  weights.   Re- 
garding this  saline  as  a  simple  beam,  the  Sernoulli-Huler 
law  can  be  stated  as 


Ll(x)  =   31  (21) 


where, 


B. 


".:(::)  is  the  bending  moment;  3  is  the  Young's   modulus 

for  the  material  of  the  beam;  I  is  the  geometric  moment  of 

inertia  which  depends  on  the  dimensions  of  the  beam  and  E(x) 

is  the  radius  of  curvature  for  the  curve  formed.   Assuming 

small  deflection 

R(x)  =   1  (22) 

Y^Tx") 

where,  y(x)  denotes  the  beam  elastica. 

In  analogy,  the  definition  of  a  mathematical  spline  can 
be  stated  in  the  words  of  Ahlberg,  Nelson  and  Walsh  (7)  as 
follows: 

"The  mathematical  saline  is  the  result  of  replacing  the 
draf tman' s  spline  by  its  elastica  and  approximating  the 
latter  by  a  piecevri.se  cubic  (normally  a  different  cubic 
between  each  r-air  of  adjacent  ducks)  wdth  certain  dis- 
countinuities  of  derivatives  permitted  at  the  junction 
points  (the  ducks  where  two  cubic  join)". 
A  mathematical  spline  is  continuous  and  has  a  continuous 
first  derivative  as  well  as  a  continuous  second  derivative. 
Thus,  s-nline  interpolating  functions  are  a  class  of  piecewise 
interpolating  polynomial  functions  satisfying  certain  con- 
tinuity TjroTserities  at  the  interpolating  points.   The  idea 
of  spline  a-oT>roximation  was  first  pointed  out  in  1946  by 
Schoenberg  (3).   In  1949,  Sard  (9)  generalised  the  classical 
abroach  by  means  of  searching  the  best  approximating  function 


?. 


of  order  1,1,  where  1  £  M  <  IT  such  that  equation  (1)  is  exact 

for  polynomials  of  degree  (II— X )  or  less.  He  then  fixed  the 

(N-I")  degrees  of  freedom  in  determining  the  coefficients  by 

effectively  requiring  that 

U*i(x)'~  dx  =  minimum  (23) 

a 

where,  U""(x)  is  the  II  ""  derivative  of  the  approximating 

function.   Schoenberg  (8)  proved  among  other  properties  that 

the  spline  interpolation  formula  is  the  "Best"  interpolation 

formula  in  the  sense  of  Sard. 

The  simplest  l:ind  of  spline  function  for  interpolation 

is  the  cubic  spline.   It  has  already  proved  to  be  of  great 

use  in  approximation  theory  and  system  identification.   Its 

properties  are  discussed  in  (7).   The  concept  is  extended  to 

curves  that  are  composed  of  segments  of  polynomials  curves 

of  an  arbitrary  degree  and  the  spline  function  so  fitted  are 

hnov.n  as  polynomial  splines.   The  theory  has  also  been  applied 

to  approximation  in  two  dimensions.   For  the  purpose  of  present 

discussion  we  will  only  consider  the  spline  approximation  as 

applied  to  differential  quadrature  method. 

1.4.2  Natural  Splines 

This  is  a  subclass  of  splines  used  by  Bellman,  Hashef 
and  Vasudevan  (3)  for  the  evaluation  of  weighting  coefficients 
in  equation  (4).   Natural  spline  of  degree  (211-1)  is  defined 
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by  the  following  condition. 

i)  S  is  a  polynomial  of  decree  2!.:-l  in  each  interval  (XT»xi+i)» 

I  s  i,  £,  .  .  .  .  ,  i\— X . 

ii)  Z   is  a  polynomial  of  decree  I.I-1  outside  the  region  (a,b). 

This  is  the  boundary  condition  requirement, 

( 2M  2 ") 
iii)  £,2f, ,S^~"   ;   are  continuous  at 


n  »"2' '"K 


iv)  S(xi)=u(xi)  for  i=l,  2,  .  .  .  .  ,K. 


An  ecuivalent  definition  in  the  sense  of  Schoenberg  is 
given  by 
i)  f°  3**(x)  ~d::     exists  and  is  minimised  subject  to  the 

following  two  conditions. 

ii)  SjS1, ,3^-'-'        are  continuous  in  (a,b) 

iii)  S(xi)=u(xi),  i=l,  2, IT 

Both  the  above  definitions  of  saline  uniquely  specify 
function  3(::)  on  the  interval  (a,b).   As  can  be  seen  these 
functions  have  strong  convergence  properties  and  for  this 
reason,  these  can  be  effectively  employed  in  the  approximate 
•crocess  of  interpolation,  integration  and  differentiation. 

1.4.3  Cardinal  Spline 

Cardinal  srline  is  also  called  "fundamental  spline"  and 
is  a  saline  for  which  exactly  one  defining  value  is  one  and 
all  others  ar3  sero  i.e.  cardinal  spline  is  a  natural  spline 
function  with  the  interpolating  conditions, 
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G,(::.)  =&i.  ,    i  =  1,  2,  .  .  .  .  ,K   (24) 

consider  cubic  spline  and  1,1=2,  -then  the  natural  spline  function 

S(x)  which  interpolates  u(::)  can  be  expressed  as 

N 
S(x)  =  X   U(x.)  C.  (::)  (25) 

i=l     "   x 

The  arbitrary  function  U(x)  can  be  approximated  as 

U(x)  =2   U(x.)  C.(x)  (26) 

i=l 

How,  by  operating  with  L  on  the  function  S(::)  of  equation 

(25),  we  get 

N 
L  S(x)  =  1  X   U(x. )   C.(x)  (27) 

i=l     x 

Since  U(x4  )  is  given  for  a  particular  value  of  xs ,  the 

operator  L  can  be  taken  inside  the  summation  and  therefore, 

ec"(27)  can  be  written  as 

w 

L  3(x)  =  X   U(x.  )1  C.  (::) 

i=l     x 

=  X   (I'C.tx))  U(x.)  (28) 

i=l 

comparing  equation  (23)  with  equation  (1), 

k±(::)   =  L  C±(x)  (23) 

The  procedure  can  be  applied  to  determine  the  weighting 
coefficients  of  differential  quadrature  method  as  follows: 

1.4.4  Computational  Scheme  for  '.veighting  Coefficients  a^. 's 

The  complete  set  of  cardinal  splines  which  span  the 


whole  space  depends  on  the  boundary  conditions  of  the  problem, 
Let  us  assume  that  the  boundary  conditions  of  the  problem 
under  consideration  are  hnown  in  the  form 

3u(x) 


b  x 


:-:=::.;   =  b-, 


(30) 


•aU(::) 


X—  X»T  ~      0^ 


(3D 


"here  the  function  U(x)  is  approximated  at  the  points  x, 
,s„  using  the  approximation 


»   •""">  » 


N 


7bU(::-i)   =  X   a,  .  U(x .),   i  =  1,  2,  .  .  . 


o  :: 


i=l 


13    3 


Defining  a  vector  Y  as 

Y  =  (U. • ,  U-,,  U?,  .  .  .  .  ,  U„.  UN' ) 


-' » 


where , 


Ui'  =  U(::) 


"1 


U'  =  U(x) 


(32) 
(33) 

(34) 
(35) 


Then,  the  (11+2)  cardinal  splines  which  completely  span 
the  st) ace  are  .given  as 


Ci(x.)  .6jLJ  ; 


ov  2 


C7,,,(x  )  =  0  ; 


Ci,(x1)  =  G2'(xN)  =  0 
V(x.)  =6id 


for 


^   —   X  ,  i-  f       •   .   •   .   ,-* 


0  =  1,  d ,  .  .  .  .  ,  x. 


(36) 


13. 


where  g.  .  is  the  hronoher  delta  defined 

for  i=j 


S-  •=  1 
13 


(37) 


o  txier  vase 

The  spline  interpolating  function  here  is 
3(::)  =  1   U(x.  )  C  (x)  +  U  '  C  (x)  +  UP'  C„ A,  (::) 


iv+1 


(33) 


Representing  the  "+2   elements   of  Y  as  Y. 's,    i= 


>  -  *  -  > 


.  .  .  .  ,K+1  raid  considering  the  I*+2  cardinal  snlines  as  C.(::), 
i=0,  1,  .  .  .  .  I'-rl  equation  (33)  can  be  written  as 


K+l 
S(s)  =  X   ?i  M*) 

1=1 

ITow,    operating  S(x)   by  1,    we   get 

TJ+1 

:(:■:)    =    J     [LC.(x)]      y 


llOl 


(39) 


(40) 


l=u 


comparison  of  equation   (40)    with  equation   (4)    gives 
a..  =     [LC^x)]  i,i=0,    1,    2,    ...    .    ,N+1  (41) 

J 

Similarly,  if  the  boundary  conditions  are  defined  as 


b*  U(x) 
b? 

2  ,  v 

a  u(x) 


-    JO 


1 


x=x. 


(42) 


(43) 


X — <*Wr   —   —5^J 


the  complete  set  of  cardinal  splines  over  the  whole  space 
can  be  defined  as 

Ci(xl:}  "  Sij  5        C.«(x.)  =  C^Cxjj)  =  0 
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00(*,)  =  0   ;  C0«(x.)  =5,, 


3 
CK+1(x.)  =  0  ;  *W*j>  =  5NJ 

i  =  1,  2,  .  .  .  .  ,K 
3  =  1,  2,  .  .  .  .  ,!!    (44) 
The  weighting  coefficients  can  be  determined  in  a  similar 

way  as  in  the  case  of  first  order  derivatives. 

1.4.5  Construction  of  cubic  spline  using  raw  data 

A  cubic  saline  in  the  interval  (0,1)  is  constructed  by 

considering  the  N+l  nesh  points  Xq,x-,  , ,x^  in  the  inter- 

val  (0,1)  such  that 

0=xQ  <  xx  <  ::2  <  <  x^b  (45) 

with  each  value  of  x,  is  associated  a  data  point  yi  which  is 
assumed  to  be  given  or  known.   Then  the  spline  S(x)  with 
respect  to  interval  0  to  b  should  have  the  following  properties 
i)  3  (x)  is  continuous,  together  with  its  first  and  second 

derivatives  on  (0  <  x  <  b) 
ii)  3  (::)  is  cubic  polynomial  within  each  subinterval 

x.  -,  <  x  <  x1?   i=l,  2,  .  .  .  .  ,H 
Analytically,  the  saline  can  be  represented  over  each 
sub-interval  (*4  -j  >-:5_)  -s 

S*(V)   -  -  ^  -  x>  +  p  (;:  -  ::i-l) 

"  x     -*i       "    -i      (46) 

Integrating  equation  (45)  twice  and  evaluating  the 


1  ^ 


constants  of  into -ration,  v;e  obtain 


on.  on,  -g- 

o  i 


::,hr    %  (;c-x.  x) 


(7,  -  -^ ) 


1 


(47) 


Prom  equation  (46)  and  (47),  it  can  be  seen  that  S  (:■:)  and 
S"  (::)  are  continuous  on  (0,b).   The  continuity  of  S  (::)  at 


~] 

recuires 

i 

"'"i-1 

+ 

ai+ai*l 

*'"i 

+ 

hi+l  „ 

"""i+1 

= 

V       _1T' 

"  i+1  i 

- 

•"" .  —    ~r .     .. 
°  1      1—1 

hi+l 

"*i 

(48) 


1.4.5.1  Corrcutational  Scheme 


Prom  equation  (47),  it  is  clear  that  calculation  of 
saline  S(x)  recuires  the  calculation  of  moments  M.(x),  i=0, 
1,  2,  .  .  .  .  ,11,   This  is  carried  out  as  follov.-s. 


btet- 


1.   Define  and  Calculate  the  folio  "ins  parameters, 


n        _ 


Nxl 


i   n.  +  n.  -, 
i    i+l 


ai  "  x   ui 


[("i+1  -  yi)/hi+l  "  (yi  "  yi-l)/hi^ 


o 

1  h.  +  h.  , 

i    i+l 


i  =  0,  1,  2,  .  .  .  .  ,N  (43) 

Step  2.   Initialise  Q0=0,Uc=0  (50) 


.0. 


il  1-       ll— i. 


then  calculate   Qn,  =  -  "k/Pv  ,  1c  =  1,  2, 


(51) 


n  =  L ii_^ ,   >  =  1,  2,  .  .  .  .  ,!: 


■D 


(52) 


Ste^  4.   The  moments  M-  can  then  be  calculated  using  the 
follov.dng  eouations: 

U.TT   —   o.T.T   V>?»   -1 

ri      i>   ;.<  — -L 


ai    1;— -1 

1!k  =  Qh  ::v+l  +  Uv  •   k  =  N-l,  .  .  .  .  2f  1 

,-  =  (do  -  °o  I!D  (53) 

o        ^ 

Thus,  the  st)line  can  be  determined  using  equation  (47). 
Theoretical  discussion  for  the  above  algorithm  is  given  in 
(7)  and  (10). 

1.5  Concluding  Remarks 

Differential  Quadrature  method  could  prove  to  be  a 
very  useful  tool  for  the  solution  of  non-linear  partial  dif- 
ferential eouations  since  it  provides  approximation  to  mar- 
tial derivatives  in  terns  of  functional  values.  Moreover, 
it  provides  a  simpler  way  of  conversion  of  partial  differen- 
tial eouations  to  ordinary  differential  equations.   In  most 
of  the  cases  boundary  conditions  are  reducible  to  simple 
algebraic  ecuations.   The  method  is  applied  to  the  solution 
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of  nacked  "bed  reactor  -oroblems  as  well  as  to  a  simple  case 
of  moving  boundary  value  problem,  The  results  are  discussed 
in  the  following  chapters.   For  the  evaluation  of  weighting 
coefficients,  the  method  of  splines  approximation  was  found 
to  be  too  involved  in  mathematics  and  could  not  be  experimented 
with  in  the  present  study.   Therefore,  the  weighting  coeffi- 
cient for  the  differential  quadrature  approximation  are  e- 
valuated  using  the  analogy  with  the  classical  quadrature  case. 


HAPT3H  II 


-J  1.2  2  — >-i-i~J.!.i  i.  J-ii-U      «JU.-_l.'ii..-.iUIi J     --ii'.— .'      _u i      ijlli^jii:...ijj 

pts  a  rtmnp    •.•TTrnt?     avt/it    T"T"rTT"'"' 

2.1  Introduction 

To  illustrate  the  method  of  previous  chapter,  let  us 
consider  the  case  of  a  chemical  reactor,  tubular  in  ion  and 
which  is  filled  with  a  packing  material.   Such  a  reactor  is 
called  a  fixed  bed  reactor.   The  transient  reaction  taking 
place  inside  the  reactor  is  of  the  form 

A  +  A  >      23  (1) 

For  the  sake  of  simplicity,  we  consider  the  case  of  an 

isothermal  reactor  only  in  this  chapter.   The  equations  for 

the  adiabatic  case  •.Till  be  dealt  with  in  the  following 

chapter.   Assuming  that  the  packing  material  has  no  influence 

on  the  reaction  taking  place  inside  the  reactor  except  its 

contribution  to  the  axial  mixing.   The  transient  equations 

for  this  type  of  reactor  can  be  given  as 

2 

1  d  P>     o  p.  -  ^  ^2  _  b  v 

*  e 
where  p  is  the  partial  pressure  of  reactant  A  in  the  inter- 
stitial fluid;  s  is  the  dimensionless  reactor  length;  t  is 
the  dimensionless  time,  and  N^  and  r  are  the  peclet  group 
and  the  reactor  rate  group  respectively  defined  as  follows 

"Pe    ^T"  (3) 
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19 


-L-2L.  (4) 


JT! 


(5) 


c  =  «   ■■»  v  ~  j 

where  Dn  is  the  average  diameter  of  the  packing  particle 
v  is  the  average  interstitial  velocity,  D  and  k  are  respective- 
ly the  effective  axial  diffusion  constant  and  the  chemical 
rate  constant,  and  x  and  ©  being  the  reactor  length  variable 
and  tine  variable  respectively.   The  boundary  conditions  for 
?(z,t)  are 


>  ( rs      -\  1      B  -  r,  t  n     _  A    +  > 


o 

^^  -      (S) 

v.There  s~  is  the  total  dimensionless  length  of  the  reactor 
and  ~°     is  the  concentration  of  A  before  entering  the  reactor. 
The  initial  condition  for  p(s,t)  is 

p(s,0)  =  p°  ,   a*  t=0,  0<2<zf     (8) 

The  above  "oroblea  was  solved  by  lee  (15)  using  the 
generalised  Nevrfcon-Haphson  method.   The  sane  numerical  data 
was  used  in  the  present  context  for  the  sake  of  comparison 
of  results. 

2.2  Differential  Quadrature  approximation: 

Using  approximation  (1)  and  (2)  of  Chapter  I,  the  partial 
derivatives  of  v   can  be  represented  as 


on 


d^iV"^  )     _  v"     P     r  ft)   -Pn-^  i  —  1         P  T' 


D"  do) 

Where  1"  is  the  nuiaber  of  selected  points  Z-,Z9,  .  .  .  ,  Z,.T 
and  a.  .'s  are  the  weighting  coefficients  and  the  notation 

p.(t)  refers  to  the  function  1  value  of  p(Z,t)  at  Z=Z  .. 

J  J 

Similarly ,  for  the  second  order  derivatives  of  p 

V       TT 


s=  "  (11) 

using  (10)  and  (11),  equation  (2)  b 9 cones 

X    a   a   P,(t)  -  Z  a   F,(t)  -  r?,2(t) 
i=i       «•'   **"«   °       i=l     «   " 


s 

!:=1 


Pe 

=   ^  pi^t)   ,   i  *  1,  2, S    (12) 

dt 

If  the  points  21,Z2»  ....  ,Z,T  are  so  selected  that 

Z-=0  and  Z--=Z.  the  boundary  conditions  (7)  and  (S)  can  be 
j.       r>   o 

.*  v^    —  C  ■  t~.  — U  ^   W-    ~~  — 

Pe  =  P  .(t)  -   3.    S   a.  .P.(t)   ,  S=Z1=0,  t>0 


>TDo    3=1 


(13) 


N 

£  a,T.  P.(t)  =  0,    Z=Z7.=Zf,   t  >  0     (14) 


.1=1 


and,  the  initial  condition  is 

pi(o)  =  ?i°   ,    i  =  1,  2,  .  .  .  .  ,N   (15) 
Squation  (12)  represents  a  set  of  K  ordinary  differen- 
tial equations  in  IT  variables  P-L(t),  P2(2),  ....  ,?T.f(t) 
subject  to  the  boundary  conditions  (15)  and  (14)  and  the 
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initial  condition  given  by  equation  (15).   Satiations  (13)  and 
(14)  are  ordinary  algebraic  equations  in  IT  variables  and 
therefore,  can  be  solbed  simultaneously  for  two  variables  say, 
P-,(t)  and  ??(t)  in  terms  of  other  (N-2)  variables  i.e.  using 
(13)  and  (14),  we  can  obtain 

?n(t)  =  i\(t,  P.,  P.,  ....  ,?r)      (16) 

?c(t)  =  fc(t,  P.,  PA,  ...  .  ,P«)      (17) 

Substituting  equations  (16)  and  (17)  into  equation  (12), 
we  can  obtain  a  set  of  (N-2)  ordinary  differential  equations 
in  (N-2)  variables 
d?,(t) 


u. 


-      ±±K   C,  £y      -4,   •   •   •   •   J-jW 


*  i  =  3,  4,  .  .  .  .  ,N    (13) 

Ecuations  (1?)  can  be  solved  with  the  given  initial  con- 
ditions (15).   Fourth  order  Range-Nutta  method  was  used  for 
this  -ournose  in  the  present  analysis. 

2.3  Computational  Procedure 

Step  1.   Select  N  appropriate  points  between  0  and  Zf  repre- 
senting the  values  of  S-,  ,  3p»  »  •  •  •  r^T  such  that  '^^=0, 

zN=sf. 

Ste^   2.      Calculate   the   weighting  coefficients   ^.'s  by  solving 
the   set  of  algebraic   equations 


1       ZT1-   a..   =    (l:-i)    zrfc    ,      l:  =  1,    -,....    , 


3=1 


v_i  .  \      1:-2 

'ij 


9  _  JT 


i   =   1,    2,    .    .    . 


>  • 


(19) 
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Step  3.   To  start,  let  t  =  0.   Initialise  P, ,  ?0,  ....  ,? 

using  equations  (15).   Let  k  =  1. 

Step  4.   Set  tT,=tv-l+At.   Solve  equations  (13)  and  (14)  for 
Ml)  and  ?0(t)  using  the  values  of  ?,(t),  PA(t),  .  .  .  .  , 

?T,T  (  t  )  . 

Step   5.      Calculate   the   co-efficients  for  Range-Eutta  method  as 

nil     =     M^V'     P^V'      ....      .W) 

ni2   =  ^^V   +  ^31'   F4(V   +  *n41'    '    '    »   W+Ha5 
21.,=  f.  (P^t-.J+ihu,,,  ?A(t1_)+-|-nAp,    .    .    fPwCtv^+^nwo) 

m.,=  f.  (?.(tnJ+n,,,   P/(t.  )+mA,f    .    .    .    ,PN(tk)+s.n) 

-    —    jj    1 1     •     •     .     •     ,i<  (20) 

Ste-  6.      p   (t,J=P,  (tlr  ,)+l    (n.1+2n,..+2ri.    +n,J, 

o 

i   =    }.    4, 


-> »     . »    •    •    •    •     >■ 


(21) 


Step.   7.      Set  !:=>+!.    and  repeat   step  4   through  7   til   a  steady- 
state   is  reached. 

2 . 4  Nume r i  c al  Re  sul t  s 

The  following  numerical  data  was  used  for  the  problem 

Pe  =  0.07  ;  Zf  =  48 

Kpe  =  20  ;  tf  =  30 

r  =  1.0  ;  P.(0)  =  v±°   =  o 
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t  =  0.2  |  i  =  1 


»  -  > 


.  .       „xa  ~tn 
Three  different  experiments  were  conducted  using  7  »  9  » 


pyi  **( 


11  ^  order  differential  quadrature  approximation.   The 
assumed  points  for  the  three  cases  are  as  follows: 

a)N=7;  C1=0,  Z2=5,  Z3«10,  Z4=20,  35=20,  S6=40,  3?=43 

b)N=9;  Zx=0,  Z2=2,  S3=5,  S4=6,  2^=10,  S6=20,  2?=30, 

Zo=40,  SQ=43 

c)H=ll;  3-,=C,  Zio=2,  ^y^*  Z.s=o,  u^=of  ig-xu,  ^-_>, 

i-g-^U,  £lg-j\J,       ^10-4Uf   -H-^° 

The  re  stilts  from  the  experiments  are  shown  in  tables  2.1, 
2.2,2.3  for  cases  (a),  (b)  and  (c)  respectively.   Pig-ares  2.1, 
2.2,  and  2.3  show  the  respective  plots  of  partial  pressure  in 
the  reactor  versus  the  dimensionless  reactor  length.   The  re- 
sults of  the  finite  difference  method  (20)  arc  plotted  in 


f i  sure  2.4, 


2.5  Discussion 

As  may  be  seen,  the  results  from  the  differential  ouad- 
rature  method  are  quite  encouraging.   The  results  from  the 
seventh  order  differential  quadrature  approximation  differ 
slightly  from  those  of  finite  difference  method  but  those 
from  differential  cuadrature  method  of  order  9  very  closely 
agree  with  the  results  of  the  finite  difference  method.  How- 
ever, finite  difference  method  requires  functional  evaluation 


24. 


for  about  480  different  points.   Thus,  besides  the  simplicity 

of  differential  quadrature  method,  major  advantage  is  obtained 
in  con-nutation  tine  which  is  much  less  than  that  required  by 
finite  difference  method.   Only  about  half  a  ninute  of  con- 
putation  tine  is  required  by  9°  *  order  differential  quadrature. 

4-V, 

The  result  from  the  11   order  nethod  are  similar  to  those 

■Hi 

from  9 u"  order  but  required  more  computation  time,   iloreover, 

instability  problems  arises  with  higher  order  methods.   Thus, 
it  can  be  concluded  that  differential  quadrature  of  order  9 
is  good  enough  for  the  solution  of  the  model  of  Isothermal 
Reactor  model  vvith  axial  miming. 


25, 


Table  2.1 

Isothermal  Reactor  with  azial  mixing  "by  7™  order  Differential 

Quadrature 


T  P(t,0)  P(t,5)  P(t,10)  P(t,20)  P(t,30)  P(t,40)  P(t,48) 

00  0  0  0  0  0  0 

5  0.0651  0.0280  C.0094  -0.0003  0.0001  -0.0000  0.0000 

10  0.0679  0.0451  0.0241  0.0024  -0.0003  0.0001  -0.0000 

20  0.0687  0.0568  0.0447  0.0197  0.0034  -0.0003  0.0050 

30  0.0675  0.0502  0.0420  0.0308  0.0166  0.0039  -0.0000 

40  0.0625  0.0292  0.0319  0.0315  0.0222  0.0158  0.0030 

50  0.0488  -0.0401  -0.0149  0.0283  0.0208  0.0219  0.0095 

60  -0.0026  -0.3410  -0.2712  -0.0229  0.0196  0.0199  0.0142 
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Isothermal  Reactor  with  Axial  Mixing  by  7 
order  Differential  Quadrature 
Z, =(0,5, 10,20,30,40,48) 
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Isothermal  Reactor  with  Axial  Mixing  by  9 

order  Differential  Quadrature 
Z, =(0,2, 5, 6, 10, 20, 30, 40, 48) 
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Isothermal  Keactor  with  Axial  Mixing  Dy  11 
order  Differential  Quadrature 
Z±=(0, 2, 5, 6, 8, 10, 15, 20, 30, 40, 48) 
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Isothermal  Reactor  with  Axial  Mixing  by 
Finite-Difference  Method 
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CHAPT3R  III 

DIFFERENTIAL  QUADRATURE  AND  THE  ADIABATIC 
RE ACTOR  WITH  AXIAL  MIXING 

3.1  Introduction 

A  more  general  case  of  pached  bed  reactor  is  the  one 
where  enerrw  balance  is  also  the  criterion  in  addition  to  mass 
balance.   Such  is  the  case  of  an  adiabatic  reactor  with  axial 
mixing.   The  differential  quadrature  method  of  Chapter  I 
■nroved  to  be  an  effective  tool  to  solve  this  problem  for 
steady  state  lihe  in  the  case  of  an  isothermal  reactor.   As- 
suming the  same  chemical  reaction  and  the  same  role  of  the 
packing  material  as  in  the  case  of  an  isothermal  reactor,  the 
equations  representing  the  dynamics  of  an  adiabatic  reactor 
with  axial  mixing  can  be  given  as  follows 

1  z2v   _  a?  _  r  v2   exu  (-3/RT)  =  3?   (1) 

1__  32T  _  3T  +  0   r0  T)2   exT>(-E/HT)   =  aT      (2) 
NPe  I?       8S  8t 

v.he  re , 

r0=DP  1:0  (3) 

v 

Q   =  K_  (4) 

cf  pf 

kn  frequency  factor  constant  in  the  Arrhenious  equation. 
AH  Heat  of  the  reaction. 
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c^  Specific  heat  of  the  reaction  mixture. 
P.P  Density  of  the  resection  mixture. 

Other  variables  have  the  sane  meaning  as  in  the  case  of 
the  isothermal  reactor.   [The  mas  axial  diffusion  co-efficient 
is  assumed  to  be  ecual  to  the  thermal  axial  diffusion  con- 
stant.  The  boundary  conditions  for  equation  (1)  and  equation 
(2)  are  respectively, 

P  =  ?(0,t)  -  _L_     a?  ,  at  3=0,  t >  0    (5) 
x'Pe  a" 


S?_  =  0    at  Z=Z~,  t>  0  (6) 

32 


T  =  T(0,t)  -  1    21  ,  at  Z=0,  t  >0   (7) 
1JPe 


a; 


3T  =  0    at  Z=Zf,  t>  0  (8) 

as 

where.  T  represents  the  temperature  of  the  reaction  mixture 
before  entering  the  reactor.  The  initial  conditions  for  (1) 
and  (2)  are 

P(Z,0)   =  P°(Z),    at   0<S<Zf,    t>0  (9) 

T(Z,0)    =   T°(Z),    at   0<Z<Zf,    t>0  (10) 

The  solution  of  the  above  problem  was  obtained  by  Liu 
and  Amundson  (15)  using  finite  difference  method  and  by  Lee  (15) 
using  the  generalised  Uevrton-Raphson  method  for  the  same  nu- 
merical data.   The  results  are  shovn  in  figures  3.13  and  3.14. 


■> 


4. 


itt't  rttiyiT 


The  solution  of  this  nroblem  was  obtained  using;  differential 
Quadrature  method  follov.-ing  the  procedure  given  in  the  next 
section. 

3.2  Problem  Formulation 

By  using  equations  (4)  and  (9)  of  Chapter  I,  the  partial 
differential  equations  for  the  adiabatic  reactor  (equations 
(1)  and  (2),  (5)-(lC))  can  be  reduced  to  a  set  of  21:  ordin; 
differential  equation.   Defining, 

?(2i,t)  =  ?i(t)   ,   i  =  1,  2,  3,  .  .  .  .  ,K 

?(Si,t)  =  T±(t)   ,   i  =  1,  2,  3,  .  .  .  .  ,N 

Equations  (1)  and  (2)  can  be  represented  as 

L?e  k=l  3=1  3=1 

r  (?, (t))2  exT3(-S  )  =  ?i(t)   ,  i=l,  2,  .  .  .  ,N 

01  RT     — ra — 

(11) 


"?e  h=l  3=1  3=1 


,  9 
•a  i" 


(P.(t))"  exiD(-E   )  =  dTi(t)   ,  i  =  1,  2,  .  .  .  ,1! 

0   x  ^T     

at  (12) 

issuminq  Zn =0   and  ZTT=Z.r,    the  boundary  conditions  of  the 
problem  can  be   represented  as 

Pe-PiC*)   "     1        1       ai3P3(t))'    aUs2r0t>° 
^e      j=l 

(13) 


15 


4   a,..  P.(t)  =  0,    at  2=Z,j=Z.p,  t  >  0   (14) 
3=1 

T  =  T-(t)  -  1   2   £ux  ?.(t)),  at  Z=Z,=0,  t  > 

x      m —  2-    - u      3  - 

"?e  j=l 

(15) 
35   a,,,  T,(t)  =  0,   at  Z=Z,T=Z^,  t  >  0    (16) 


j=l 


which  gives  a  set  of  four  algebraic  eo_uations.   Similarly, 
the  initial  conditions  (9)  and  (10)  can  be  represented  as 
?i(3,0)  =  ?i°(Z),   i  =  1,  2,  .  .  .  .  ,11,  t  >  0 

(17) 
^(2,0)  =  Ti°(S),   i  =  1,  2,  .  .  .  .  ,N,  t  >  0 

(18) 
She  set  of  equations  (11)  and  (12)  can  be  solved  subject 
to  boundary  conditions  (13)-(17)  and  initial  conditions  (18) 
and  (19)  by  following  the  computational  procedure  as  given 
belov;. 

3.3  Computational  Procedure 

Step  1.   Select  IT  and  Z^,i=l,  2,  .  .  .  .  ,N. 

Ste-n  2.   Determine  the  weighting  co-efficients  a.  -,i=l,  2, 
....  ,IT  and  3=1,  2,  ...  .  ,N  by  using  the  selected 
values  of  Z.  and  solving  the  set  of  following  algebraic 


eouations 

=  (k-D  z. 


z         A.     a..  =  (L-l)  Z«     ,  L-  1,  2,  .  .  .  ,H 


j=l  i  =  1,  2,      .  ,N 

(19) 
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Step  3.  Solve  equations  (13)-(16)  to  obtain 

?1(t)  =  i-l  (?3,  P4 ,?N)      .  (20) 

?2(t)  =  f2(P3,  ?4 ,PN)         (21) 

E1(t)  =  P1(T3,  T4,  ....  ,Tn)         (22) 

T2(t)  =  P2(53,  T4,  .  .  .  .  ,Tn)         (23) 

Step  4.   Substitute  equations  (20)-(23)  in  (11)  and  (12)  to 
get  a  set  of  (2IT-4)  ordinary  differential  equations 

a  M*)    -  f  '  rp     -d  .p..  t.). 

C.  u 

i  =  3,  4,  .  .  .  .  ,H    (24) 

dT. (t)    —  T?  '   f»P    T  «p-p^ 

1  dt        1    3   4  ..   i 

i  =  3,  4,  .  .  .  .  ,»    (25) 

Step  5.   Solve  equations  (24)  and  (25)  rath,  initial  conditions 
given  by  equations  (17)  and  (18). 

Step  6.   Calculate  P1,?2,T-L,T2,  t  >  0  by  using  the  results 
from  step  5  and  equations  (20)-(23). 

Step  7.  Repeat  steps  3-6  till  a  steady  state  is  reached. 

3.4  Numerical  Results 

The  numerical   data  used  was  taken  from  that  used  by 
Lee    (15)   and  also  by  Liu  and  Amundson  (15)   and  is  given  below: 
Pe  =   0.07   atm.         ;        3/R  =   22,000°R 


37. 


T  =  1250°?.     ;    Q  =  0.1  ::  104  °R/atm. 


e 
!L  =2.0       ;    r  =  0.5   x  10S  (ata)"1 


=  0.2        ;     Zx-  =  43 


.lb    l/=W 


P°(Z)  =  0       ;    T°(Z)  =  1270°P. 
A  number  of  experiments  were  done  in  order  to  analyse 
the  behavior  of  these  equations.   The  order  of  differential 
quadrature  approximation  used  and  the  selected  points  in 
each  case  is  given  below: 

a)  K=7;      Z-^0,    Z2=5,    Z3=10,    Z4=20,    2^=30,    Zg=40, 

Z7=48 

b)  N=9;      S1=0,    Z2=2,    Z3=5,    S4=6,    Z5=10,    Z6=20, 

Z7=30,    Zfi=40,    ZQ=48 

c)  N=ll;    Z-^0,    Z2=2,    Z3=5,    Z4=6,    Z5=3,    Z6=10, 

^rj=±Sf         £g-<-U|        iQ-jU,         £^Q—ifvJ,         ^T_2.  —  *W 

d)  n=9;     2^=0,  z2=o.5,  z3=i.o,  z4=5.o,  z5=io, 

Z6=20,    Z7=30,    Z8=40,    Zg=48 

e)  N=9;      31=0,    Z2=0.1,    Z3=1.0,    Z^ag.0, 

7        in        n     —OO        7     —  "5  0        7-—  ^  fi        7     =ZLfi 

■r\    w   to.    7    _n      7    —  r>    -      7    — ^      ^    — 1 0      "    -l  ^      7.^=20. 

I  J     I\=lii;     _.-,=0,     j9  =  0.>f     Zi  -j  —  •,     ±i£—  ±'Jf     tie— 13|     ■''6        u» 

^=0,    og=3^,    «g=j-/f    ~io    ■",    iji]_-^-'»    "12~^W 


33. 


rr\    TJ— 1  ?•     7    —0        7,     — D    R  7.    —  1         7     —  R        7    —1  Pi        7    _ T  c 

—  —  J  *f  ^  C 

S?=20,    z8=25,  Zg=30,    S10=35,    Si:L=40, 

h)   I;=9;      Z^O,    Z2=0.5,  Z^l,    24=5,    25=15,    Sg=25, 

•7       _■)'  7       —  /  CJ  "       —fir, 

;Jr-r-J>,       iJg  — <r>,  ZJq-^-U, 

■i  ^     "— O  •         7     _n        7     —C     R  7     _9        7     —  £        7     — T  fi        7     _on 

zr- 


=30,    SP=39,    SQ=43 


The  results  iron  experiments  (a)-(e)  and  (i)  ere  tabu- 
lated in  tables  3.1  -  3.12  respectively  and  the  respective 
plots  are  shorn  in  fig.  3.1  -  3.12.   The  results  of  the  ap- 
proximation for  the  rest  of  the  experiments  were  found  to  be 
quite  unstable  and  in  the  first  few  iterations  only,  the 
fluctuations  were  high  enough  to  exceed  the  overflow  limits 
of  the  computer.   The  result  obtained  from  the  finite  dif- 
ference method  are  plotted  in  figures  3.13  2-nd  3.14. 

3.5  Discussion 

As  may  be  seen  the  same  distribution  of  node  points  as 
in  the  case  of  isothermal  reactor,  did  not  prove  to  be  equally 
good  here.   In  fact,  the  selection  of  node  points  plays  a  very 
important  role  in  the  case  of  adiabatic  reactor.   A  seven 

point  differential  quadrature  approximation  was  not  found 

th 

to  be  satisfactory.   The  results  obtained  by  using  7  "  order 
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differential  quadrature  method  (fig.  3.1  and  3.2)  do  not  agree 
with  those  from  the  finite  difference  method.   Therefore,  the 
number  of  points  were  increased  to  9  in  further  experimenta- 
tion.  In  experiments  (a)-(b)  and  (c),  the  initial  boundary 
slope  of  the  curve  obtained  changes  with  time,  as  a  result 

of  which  wide  fluctuations  are  produced  as  the  reactor  length 

■fcVi 

increases.   Experimentation  was  also  done  using  11   order 

tl~ 
and  12 u  "  order  differential  quadrature  method  (exp.  c,  f,  g) 

but  this  was  not  found  to  be  of  much  help,  due  to  the  increased 
round  off  and  truncation  errors  in  the  computations.   Therefore, 
the  9  ~L   order  differential  quadrature  method  was  used  and  the 
initial  points  were  kept  close  to  each  other  (cases  d,  e, 
and  i)  in  order  to  maintain  a  constant  slope.   Variations 
were  considered  based  on  the  experimental  judgement.   Steady 
state  conditions  were  best  obtained  in  experiment  (i).   The 
results  obtained  (fig.  3.11  and  fig.  3.12)  agree  very  closely 
with  those  from  the  finite  difference  method  (fig.  3.13  snd 
fig.  3.14).   The  computation  time  required  was  found  to  be 
less  than  a  minute  while  it  requires  several  minutes  of  corn- 
nutation  time  by  the  finite  difference  method.   Thus,  it  can 
be  concluded  that  with  the  proper  selection  of  node  points, 

■fcVi 

9   order  differential  quadrature  method  gives  as  accurate 
results  as  the  finite  difference  method  does. 


40. 


Table  3.1 

Transient  in  Partial  Pressure,  Adiabatic  Reactor  by  7 
order  Differential  Quadrature 


T  P(t,0)  P(t,5)  P(t,10)  P(t,20)  P(t,30)  P(t,40)  P(t,48) 

00  0  0  0  0  0  0 

5  0.06*19  0.0272  0.0093  -0.0003  0.0001  -0.0000  0.0000 

10  0.0677  0.0435  0.023^  0.0023  -0.0003  0.0001  -0.0000 

20  0.0677  0.0503  0.0382  0.0175  0.0032  -0.0002  0.0000 

30  0.0664  0.0411  0.0299  0.0221  0.0136  0.0036  0.0000 

40  0.0633  0.0292  0.0247  0.0190  0.0137  0.0124  0.0026 

50  0.0564  -0.0050  0.0028  0.0191  0.0106  0.0128  0.0070 
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Table  3-2 
Transient  in  Temperature,  Adiabatic  Reactor  by  7th  order 
Differential  Quadrature 


t 

T(t,0) 

T(t,5) 

T(t,10) 

T(t,20) 

T(t,30) 

T(t,40) 

T(t,48) 

0 

1270 

1270 

1270 

1270 

1270 

1270 

1270 

5 

1251.5 

1262.4 

1267.5 

1270.1 

1270.0 

1270.0 

1270.0 

10 

1251.1 

1260.0 

1265.3 

1269.5 

1270.0 

1270.0 

1270.0 

20 

1254.2 

1277.7 

1280.5 

1271.7 

1269.4 

1270.0 

1270.0 

30 

1254.5 

1288.5 

1302.9 

1295.3 

1276.4 

1269.6 

1270.0 

40 

1249.0 

1259-9 

1235.6 

1314.3 

1301.7 

1279-9 

1270.4 

50 

1223.5 

1130.5 

1191.6 

1300.2 

1314.5 

1310.0 

1278.9 

42. 
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Transient  in  Partial  Pressure,  Adiabatic  Keactor  with 
Axial  Mixing  by  7th  order  Differential  Quadrature 
Z±={Ot 5, 10, 20, 30,40,48) 
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53. 


Transient  in  Temperature,  Adiabatic  Keactor  with  Axial 
Mixing  by  7   order  Differential  Quadrature 
2^(0,5,10,20,30,40,48) 


m 

•p 
■p 

CO 

a 


1320  _ 


1310 


1300  - 


1290  - 


1 280  - 


1270   - 


1260   - 


1250 


10      20      30      40 

Keactor  Length,  z  »■ 

.figure  3.2 


54. 


Transient  in  Partial  Pressure,  Adiabatic  Reactor  with 
Axial  Mixing  by  9   order  Differential  Quadrature 

Z.=(u,2,5,6,10,20,30,40,48) 
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55. 


Transient  in  Temperature,  Adiabatic  Keactor  with  Axial 
Mixing  by  ytl:1  order  Differential  Quadrature 
Z^Co, 2, 5, 6, 10,20, 30, 40, 48) 
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Transient  in  Partial  Pressure,  Adiabatic  Keactor  with 
Axial  Mixing  by  11th  order  Differential  Quadrature 

Z. =(0,2, 5, 6, 8, 10, 15, 20, 30, 40, 48) 
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Transient  in  Temperature,  Adiabatic  Reactor  with  Axial 
Mixing  by  11th  order  Differential  Quadrature 
2^(0,2,5,6,8,10,15,20,30,4-0,48) 
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Transient  in  Partial  Pressure,  Adiabatic  Keactor  with 
Axial  Mixing  by  yth  order  Differential  Quadrature 
2^(0,0.1,1.5,10,20,30,40,48) 
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Transient  in  Temperature,  Adiabatic  Keactor  with  Axial 
Mixing  by  9   order  Differential  Quadrature 
Z^Co, 0.1, 1,5, 10, 20, 30, 40, 48) 
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Transient  in  Partial  Pressure,  Adiabatic  Reactor  with 
Axial  Mixing  by  9   order  Differential  Quadrature 
Z. =(0,0. 5, 1,5,10,20, 30,40,48) 
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Transient  in  Temperature,  Adiabatic  Reactor  with  Axial 
Mixing  by  9   order  Differential  Quadrature 
Z^Co, 0.5, 1.5,10,20,30,40,48) 
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Transient  in  Partial  Pressure,  Adiabatic  Reactor  with 
Axial  Mixing  by  9th  order  Differential  Quadrature 
Z^C  0,0. 8, 2, 6, 10, 20, 30, 39, 48) 
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Transient  in  Temperature,  Adiabatic  Reactor  with  Axial 

Mixing  by  9   order  Difi'erentiai  Quadrature 
Z^CO, 0.8, 2,  6, 10, 20, 30, 39, 48) 
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Transient  in  Partial  Pressure,  Adiabatic  Reactor  with 
Axial  Mixing  by  J?inite-Di±*i'erence  Method 
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Transient  in  Temperature,  Adiabatic  neactor  with  Axiai 
Mixing  Dy  Finite-Difference  Method 
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CHAPTZR  IV 
CONCLUSIONS 

Differential  quadrature  was  found  to  be  an  excellent 
tool  for  the  solution  of  non-linear  partial  differential 
equations.  Results  for  both,  the  isothermal  reactor  as  well 
as  the  adiabatic  reactor  are  encouraging  in  this  respect. 
Although,  much  experimentation  was  required  in  order  to  arrive 
at  accurate  results  but  it  is  mainly  due  to  the  values  of  the 
weighting  coefficients  used  in  the  approximation.  But  in 
practice,  these  coefficients  need  to  be  determined  only  once 
for  a  -oerticular  tyne  of  -croblems.   Therefore,  if  the  optimal 
values  are  known,  the  method  can  be  conveniently  applied. 
Besides  its  advantages  in  terms  of  savings  in  computational 
time  and  computer  storage,  the  comparison  of  the  method  with. 
finite-difference  method  shows  that  differential  quadrature 
is  much  easier  to  arj-oly  in  practice.   Differential  quadrature 
using  s-nline  approximation  may  enhance  its  chances  of  success 
and  may  enable  the  successful  solution  of  more  general  cases 
of  packed-bed  reactors  without  any  instability  problems.  Thus, 
more  research  in  this  field  may  prove  to  be  beneficial. 
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ABSTRACT 

Differential  quadrature  is  a  useful  numerical  technique 
for  solving  non-linear  partial  differential  equations.   It 
involves  approximating  the  partial  derivatives  by  a  linear 
combination  of  functional  values  and,  therefore,  provides 
an  easy  method  of  transformation  of  partial  differential 
equations  into  a  set  of  ordinary  differential  equations. 
The  technique  is  employed  for  solving  boundary  value  problems 
which  can  be  represented  by  partial  differential  equations. 

Host  other  methods  like  the  finite-difference  method 
involve  approximation  in  terms  of  functional  differences 
instead  of  functional  values  and  therefore,  require  functional 
evaluation  at  a  large  number  of  points  for  satisfactory  re- 
sults.  It  is  in  this  respect  that  differential  quadrature 
has  its  major  advantages  over  other  methods  in  terms  of  both, 
the  computer  storage  and  computational  time.   However,  the 
success  of  the  method  depends  largely  upon  the  method  of 
evaluation  of  weighting  coefficients.   Three  methods  are  con- 
sidered in  this  respect  viz. classical  quadrature  analogy, 
Legendre  polynomial  approach  and  spline  approximation. 

Differential  quadrature  is  applied  to  solving  several 
models  in  engineering  with  both  fixed  and  moving  boundary 
conditions.   A  moving  boundary  condition  is  specified  at  a 
T3oint  which  itself  varies  as  a  function  of  time.   Differential 


Quadrature  is  used  to  solve  the  isothermal  reactor  model  as 
well  as  the  adiabatic  reactor  model.   A  lot  of  computer  memory 
and  commutation  time  are  saved  by  using  this  technique. 


